Congreso SP(2)

Presentación principales resultados. Parte 2.

Author

Andrés González Santa Cruz

Published

October 14, 2025


Cargar datos y exploración

Cargar paquetes estadísticos y gestión de BBDD

Code
# invisible("Only run from Ubuntu")
# if (!(Sys.getenv("RSTUDIO_SESSION_TYPE") == "server" || file.exists("/.dockerenv"))) {
#   if(Sys.info()["sysname"]!="Windows"){
#     Sys.setenv(RETICULATE_PYTHON = "/home/fondecytacc/.pyenv/versions/3.11.5/bin/python")
#   }
# }

#clean enviroment
rm(list = ls()); gc()

time_before_dedup2<-Sys.time()

if(!require(reticulate)){install.packages("reticulate")}

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

# Busca .mamba_root/envs/py311/python.exe desde getwd() hacia padres
find_python_rel <- function(start = getwd(),
                            rel = file.path(".mamba_root","envs","py311","python.exe")) {
  cur <- normalizePath(start, winslash = "/", mustWork = FALSE)
  repeat {
    cand <- normalizePath(file.path(cur, rel), winslash = "/", mustWork = FALSE)
    if (file.exists(cand)) return(cand)
    parent <- dirname(cur)
    if (identical(parent, cur)) return(NA_character_)  # llegó a la raíz
    cur <- parent
  }
}
# --- Bootstrap reticulate con ruta relativa a getwd() ---
if(Sys.info()["sysname"]!="Windows"){
  #Sys.setenv(RETICULATE_PYTHON = "usr/bin/python3")
  reticulate::py_config()
} else {
  py <- find_python_rel()
  if (is.na(py)) {
    stop("No se encontró Python relativo a getwd() (buscando '.mamba_root/envs/py311/python.exe').\n",
         "Directorio actual: ", getwd())
  }
  # Forzar ese intérprete
  Sys.unsetenv(c("RETICULATE_CONDAENV","RETICULATE_PYTHON_FALLBACK"))
  Sys.setenv(RETICULATE_PYTHON = py)
  reticulate::use_python(py, required=T)
  reticulate::py_config()  # verificación
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

resolve_doc_dir <- function() {
  # 1) Interactivo en RStudio (mejor para desarrollo)
  if (interactive() && requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) {
    path <- tryCatch(rstudioapi::documentPath(), error = function(e) NULL)
    if (!is.null(path) && nzchar(path)) return(normalizePath(dirname(path)))
  }

  # 2) Render / knitr (Quarto)
  in_knitr <- tryCatch(knitr::current_input(), error = function(e) NULL)
  if (!is.null(in_knitr) && nzchar(in_knitr)) return(normalizePath(dirname(in_knitr)))

  # 3) Si nada funciona, devolver error claro (no fallback silencioso)
  stop("No se pudo determinar la carpeta del documento. Asegura que estás renderizando con Quarto o usando RStudio.")
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

if(grepl("ubuntu",pak::system_r_platform())){
  if (rstudioapi::isAvailable()) {
    # Code that needs RStudio
    wdpath<- dirname(dirname(rstudioapi::documentPath()))
  } else {
    # Code for non-RStudio environments (e.g., command line, server)
    # This part should use a portable method like knitr::current_input() or getwd()
    wdpath <- dirname(resolve_doc_dir())
  }
}
load(paste0(wdpath,"/congresosp.RData"))
          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  602966 32.3    1320207 70.6   835224 44.7
Vcells 1138441  8.7    8388608 64.0  1930478 14.8
python:         /home/andres/.cache/R/reticulate/uv/cache/archive-v0/9jodbrlaYfTdFRuQinKW1/bin/python3
libpython:      /home/andres/.cache/R/reticulate/uv/python/cpython-3.11.13-linux-x86_64-gnu/lib/libpython3.11.so
pythonhome:     /home/andres/.cache/R/reticulate/uv/cache/archive-v0/9jodbrlaYfTdFRuQinKW1:/home/andres/.cache/R/reticulate/uv/cache/archive-v0/9jodbrlaYfTdFRuQinKW1
virtualenv:     /home/andres/.cache/R/reticulate/uv/cache/archive-v0/9jodbrlaYfTdFRuQinKW1/bin/activate_this.py
version:        3.11.13 (main, Jul 11 2025, 22:43:55) [Clang 20.1.4 ]
numpy:          /home/andres/.cache/R/reticulate/uv/cache/archive-v0/9jodbrlaYfTdFRuQinKW1/lib/python3.11/site-packages/numpy
numpy_version:  2.3.3

NOTE: Python version was forced by py_require()
Code
#https://github.com/rstudio/renv/issues/544
#renv falls back to copying rather than symlinking, which is evidently very slow in this configuration.
renv::settings$use.cache(FALSE)

#only use explicit dependencies (in DESCRIPTION)
renv::settings$snapshot.type("implicit")

#check if rstools is installed
if(Sys.info()["sysname"]=="Windows"){
try(installr::install.Rtools(check_r_update=F))
}

check_quarto_version <- function(required = "1.7.29", comparator = c("ge","gt","le","lt","eq")) {
  comparator <- match.arg(comparator)
  current <- package_version(paste(unlist(quarto::quarto_version()), collapse = "."))
  req     <- package_version(required)

  ok <- switch(comparator,
               ge = current >= req,
               gt = current >  req,
               le = current <= req,
               lt = current <  req,
               eq = current == req)

  if (!ok) {
    stop(sprintf("Quarto version check failed: need %s %s (installed: %s).",
                 comparator, required, current), call. = FALSE)
  }
  invisible(TRUE)
}

check_quarto_version("1.7.29", "ge") 

#change repository to CL
local({
  r <- getOption("repos")
  r["CRAN"] <- "https://cran.dcc.uchile.cl/"
  options(repos=r)
})

if(!require(pacman)){install.packages("pacman");require(pacman)}

Cargando paquete requerido: pacman

Code
if(!require(pak)){install.packages("pak");require(pak)}

Cargando paquete requerido: pak

Code
pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes

No 00LOCK detected in: /home/andres/Escritorio/SER2025/renv/library/linux-ubuntu-noble/R-4.4/x86_64-pc-linux-gnu No 00LOCK detected in: /opt/R/4.4.1/lib/R/library

Code
if(Sys.info()["sysname"]=="Windows"){
if (getRversion() != "4.4.1") { stop("Requires R version 4.4.1; Actual: ", getRversion()) }
}

#check docker
check_docker_running <- function() {
  # Try running 'docker info' to check if Docker is running
  system("docker info", intern = TRUE, ignore.stderr = TRUE)
}
if(Sys.info()["sysname"]=="Windows"){
  install_docker <- function() {
    # Open the Docker Desktop download page in the browser for installation
    browseURL("https://www.docker.com/products/docker-desktop")
  }
  # Main logic
  if (inherits(try(check_docker_running(), silent = TRUE), "try-error")) {
    liftr::install_docker()
  } else {
    message("Docker is running.")
  }
}

# if(Sys.info()["sysname"]!="Windows"){
#   system("sudo apt-get update && sudo apt-get install -y libudunits2-dev libgdal-dev libgeos-dev libproj-dev")
#   system("sudo apt-get install -y libudunits2-dev libgdal-dev libgeos-dev libproj-dev")
# }

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#PACKAGES#######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

unlink("*_cache", recursive=T)

# Load or install required packages using pacman::p_load
# This handles installation, loading, and avoids duplicates
# Note: INLA requires separate installation from its dedicated repository.
pacman::p_load(
    # Package Management
    pacman,       # Streamlines loading and installing packages simultaneously
    
    # Core Data & Visualization (Part of tidyverse)
    tidyverse,    # Collection of packages for data manipulation, visualization, and more
    dplyr,        # Grammar of data manipulation (e.g., filter, mutate, group_by)
    tidyr,        # Tools for tidying data, reshaping (e.g., pivot_wider, pivot_longer)
    ggplot2,      # Grammar of graphics for powerful data plotting
    stringr,      # Consistent wrappers for common string manipulation tasks
    purrr,        # Tools for working with functions and vectors/lists (e.g., mapping)
    viridisLite,  # Provides the perceptually uniform 'viridis' color palettes
    scales,       # Graphical scales for ggplot2 (e.g., custom breaks, formatting labels for axes)
    
    # Survival & Epidemiology
    mexhaz,       # Flexible parametric hazard regression models for survival analysis
    relsurv,      # Core package for **Relative Survival Analysis** in population-based studies
    survminer,    # Extends 'survival' for high-quality Kaplan-Meier plots and survival tables
    popEpi,       # Tools for standardized mortality ratios (SMRs) and population-level survival measures
    epitools,     # Epidemiological tools for data analysis (e.g., rate ratios, exact confidence intervals)
    
    # Spatial Analysis
    sf,           # Simple Features: for handling and analyzing **vector-based spatial data** (points, lines, polygons)
    spdep,        # Spatial dependence: for analyzing spatial autocorrelation and spatial regression models
    geodata,      # Access to global spatial data (e.g., administrative boundaries, climate data)
    INLA,         # Integrated Nested Laplace Approximations (Bayesian spatial modeling) - install separately
    
    # Data Cleaning & Reporting
    janitor,      # Simple tools for examining and cleaning dirty data (e.g., cleaning column names)
    tableone,     # Creates "Table 1" summaries for descriptive statistics (baseline characteristics)
    kableExtra,   # Enhances 'knitr::kable()' for complex and professional R Markdown tables
    pander,       # A tool for rendering R objects into Pandoc markdown formats
    DT,           # Creates interactive HTML data tables (useful for detailed tables in R Markdown/Quarto)
    rio,          # Simplifies data import/export with a consistent interface ('import()')
    quarto,       # For rendering the R Markdown/Quarto poster documents (.qmd files)
    
    # Statistical & Computing Tools
    biostat3,     # Functions and datasets for biostatistics teaching and research
    coin,         # Conditional inference procedures for robust hypothesis testing (e.g., permutation tests)
    metafor,      # For heterogeneity testing and general **meta-analysis** methods
    parallel,     # Base R package for **parallel computing** (useful for speeding up bootstrap and resampling)
    cowplot,      # Streamlined plot theme and functions for arranging multiple ggplot2 plots
    grid,         # Base R package for low-level graphics (used indirectly for plot arrangement)
    ellmer,       # Commonly used for likelihood estimation or generalized mixed-effects models (LMM/GLMM)
    
    # Development & Reproducibility Tools
    devtools,     # Tools to make package development easier (often used to install from GitHub)
    installr,     # Functions for self-updating R and R packages
    liftr,        # A framework for creating and managing reproducible research projects
    pak,          # A fast, dependency-aware package installer (alternative to install.packages)
    
    install = TRUE
)

# ----------------------------------------------------------------------
# 4. BPMN from GitHub (not on CRAN, so install via devtools if missing)
# ----------------------------------------------------------------------
if (!requireNamespace("bpmn", quietly = TRUE)) {
  devtools::install_github("bergant/bpmn")
}


#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#FUNCTIONS######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#
# NO MORE DEBUGS
options(error = NULL)        # si antes tenías options(error = recover) o browser)
options(browserNLdisabled = FALSE)


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#
#NAs are replaced with "" in knitr kable
options(knitr.kable.NA = '')

pander::panderOptions('big.mark', ',')
pander::panderOptions('decimal.mark', '.')

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#

#to format rows in bold
format_cells <- function(df, rows ,cols, value = c("italics", "bold", "strikethrough")){
  # select the correct markup
  # one * for italics, two ** for bold
  map <- setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough"))
  markup <- map[value]  
  for (r in rows){
    for(c in cols){
      # Make sure values are not factors
      df[[c]] <- as.character( df[[c]])

      # Update formatting
      df[r, c] <- ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup))
    }
  }
  return(df)
}
#To produce line breaks in messages and warnings
knitr::knit_hooks$set(
   error = function(x, options) {
     paste('\n\n<div class="alert alert-danger" style="font-size: small !important;">',
           gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),
           '</div>', sep = '\n')
   },
   warning = function(x, options) {
     paste('\n\n<div class="alert alert-warning" style="font-size: small !important;">',
           gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),
           '</div>', sep = '\n')
   },
   message = function(x, options) {
     paste('<div class="message" style="font-size: small !important;">',
           gsub('##', '\n', x),
           '</div>', sep = '\n')
   }
)

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#CONFIG #######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

options(scipen=2) #display numbers rather scientific number

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

# Define the function first
#joins these values with semicolons and optionally truncates the result if it exceeds a specified width.
toString2 <- function(x, width = NULL, ...) {
    string <- paste(x, collapse = "; ")
    if (missing(width) || is.null(width) || width == 0) 
        return(string)
    if (width < 0) 
        stop("'width' must be positive")
    if (nchar(string, type = "w") > width) {
        width <- max(6, width)
        string <- paste0(substr(string, 1, width - 3), "...")
    }
    string
}
normalize_txt <- function(x) {
  x|>
    stringi::stri_trans_general("Latin-ASCII")|>
    tolower()|>
    trimws()
}

# Custom function for sampling with a seed
sample_n_with_seed <- function(data, size, seed) {
  set.seed(seed)
  dplyr::sample_n(data, size)
}

# Function to get the most frequent value 
most_frequent <- function(x) { 
  uniq_vals <- unique(x)
  freq_vals <- sapply(uniq_vals, function(val) sum(x == val))
  most_freq <- uniq_vals[which(freq_vals == max(freq_vals))]
  
  if (length(most_freq) == 1) {
    return(most_freq)
  } else {
    return(NA)
  }
}

sum_dates <- function(x){
  cbind.data.frame(
    min= as.Date(min(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01"),
    p001= as.Date(quantile(unclass(as.Date(x)), .001, na.rm=T), origin = "1970-01-01"),
    p005= as.Date(quantile(unclass(as.Date(x)), .005, na.rm=T), origin = "1970-01-01"),
    p025= as.Date(quantile(unclass(as.Date(x)), .025, na.rm=T), origin = "1970-01-01"),
    p25= as.Date(quantile(unclass(as.Date(x)), .25, na.rm=T), origin = "1970-01-01"),
    p50= as.Date(quantile(unclass(as.Date(x)), .5, na.rm=T), origin = "1970-01-01"),
    p75= as.Date(quantile(unclass(as.Date(x)), .75, na.rm=T), origin = "1970-01-01"),
    p975= as.Date(quantile(unclass(as.Date(x)), .975, na.rm=T), origin = "1970-01-01"),
    p995= as.Date(quantile(unclass(as.Date(x)), .995, na.rm=T), origin = "1970-01-01"),
    p999= as.Date(quantile(unclass(as.Date(x)), .999, na.rm=T), origin = "1970-01-01"),
    max= as.Date(max(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01")
  )
}

is_stata_ok <- function(x) {
  nchar(x) <= 32 & grepl("^[A-Za-z][A-Za-z0-9_]*$", x)
}

to_ascii_lower <- function(x) stringi::stri_trans_general(x, "Latin-ASCII")

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# ── Tidy RMS ────────────────────────────────────────────────────────

tidy_cph <- function(model) {
  if (!inherits(model, "cph")) stop("Model must be an rms::cph object")
  # run summary
  s <- summary(model)
  df_s <- as.data.frame(unclass(s))
  df_s$rown <- rownames(s)
  # find rows that are "Hazard Ratio"
  hr_rows <- trimws(df_s$rown) == "Hazard Ratio"
  # build tidy tibble
  tibble::tibble(
    term     = trimws(rownames(s)[which(hr_rows) - 1L]), # previous row = variable name
    estimate = df_s$Effect[hr_rows],                     # Hazard Ratio
    conf.low = df_s$`Lower 0.95`[hr_rows],
    conf.high= df_s$`Upper 0.95`[hr_rows],
    p.value  = df_s$P[hr_rows]
  )
}

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# ── SMR Phi ────────────────────────────────────────────────────────


sir_ci_phi_improved <- function(sir_obj, phi, conf.level = 0.95) {
  #Método log-normal, the best, dont overestimate, or subestimate variance
  # extract totals
  total_obs <- sir_obj$observed
  total_exp <- sir_obj$expected
  
  # Calculate SEs
  theta <- total_obs / total_exp
    # Normal approximation, n>20  
  # Corrected SEs (McCullagh & Nelder, 1989)
  # “For ratios of Poisson means (such as SIR or CMR), the appropriate approach is to use multinomial or binomial models conditioned on the total observed.”
  # Breslow NE, Day NE. Statistical Methods in Cancer Research, Vol. II (IARC, 1987), §2.2. – Derives the same SE formula and recommends inflating by φ in the presence of heterogeneity.
  z <- qnorm(1 - (1 - conf.level)/2)

  se_log <- sqrt(phi / total_obs)  # Valid formula
  
  # ICs
  lci <- theta * exp(-z * se_log)
  uci <- theta * exp(z * se_log)
  
  data.frame(
    SIR = theta,
    CI_low = lci,
    CI_high = uci,
    phi_used = phi
  )
}

smr_print <- function(sir_out, name = NULL) {
    # Validate input
    if (!is.list(sir_out)) stop("sir_out must be a list or list-like object.")
    required <- c("observed", "pyrs", "expected", "sir")
    missing_req <- setdiff(required, names(sir_out))
    if (length(missing_req)) stop("sir_out is missing required elements: ", paste(missing_req, collapse = ", "))
    
    # Pull components (allow vectors)
    obs <- as.numeric(sir_out$observed)
    pyrs <- as.numeric(sir_out$pyrs)
    expc <- as.numeric(sir_out$expected)
    sir   <- as.numeric(sir_out$sir)
    sir_lo <- if ("sir.lo" %in% names(sir_out)) as.numeric(sir_out$sir.lo) else rep(NA_real_, length(sir))
    sir_hi <- if ("sir.hi" %in% names(sir_out)) as.numeric(sir_out$sir.hi) else rep(NA_real_, length(sir))
    ear    <- if ("EAR" %in% names(sir_out)) as.numeric(sir_out$EAR) else rep(NA_real_, length(sir))
    
    n <- max(length(obs), length(pyrs), length(expc), length(sir), length(sir_lo), length(sir_hi), length(ear))
    rep_len <- function(x) if (length(x) == n) x else rep(x, length.out = n)
    
    obs   <- rep_len(obs)
    pyrs  <- rep_len(pyrs)
    expc  <- rep_len(expc)
    sir   <- rep_len(sir)
    sir_lo <- rep_len(sir_lo)
    sir_hi <- rep_len(sir_hi)
    ear   <- rep_len(ear)
    name  <- if (is.null(name)) rep(NA_character_, n) else rep_len(as.character(name))
    
    # Format SMR: "x.xx (y.yy–z.zz)" or "x.xx" if bounds missing
    smr_fmt <- vapply(seq_len(n), function(i) {
        if (!is.na(sir_lo[i]) && !is.na(sir_hi[i])) {
            sprintf("%.2f (%.2f–%.2f)", sir[i], sir_lo[i], sir_hi[i])
        } else {
            sprintf("%.2f", sir[i])
        }
    }, FUN.VALUE = character(1), USE.NAMES = FALSE)
    
    # Format EAR
    ear_fmt <- ifelse(is.na(ear), NA_character_, sprintf("%.2f", ear))
    
    # Build result
    out <- data.frame(
        total    = name,
        observed = round(obs, 0),
        pyrs     = round(pyrs, 0),
        expected = round(expc, 0),
        SMR      = smr_fmt,
        EAR      = ear_fmt,
        stringsAsFactors = FALSE,
        row.names = NULL
    )
    
    out
}


# ── Plots ────────────────────────────────────────────────────────

theme_sjPlot_manual <- function() {
  theme(
    plot.title = element_text(size = 16, face = "bold"),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    panel.background = element_rect(fill = "white"),
    panel.grid.major = element_line(color = "gray80"),
    panel.grid.minor = element_line(color = "gray90")
  )
}


# --- Paleta base continua (tu cont_palette) ---
andal_pal <- function(reverse = FALSE) {
  cols <- c(
    "#052709", # dark teal
    "#5FA9B3", # teal accent
    "#A17C6C", # warm mid
    "#D2A899", # warm light
    "#906F00", # earth gold
    "#527700"  # olive
  )
  if (reverse) cols <- rev(cols)
  grDevices::colorRampPalette(cols, interpolate = "spline")
}

# Continua (equivalente a scale_fill_viridis_c)
scale_fill_andal_c <- function(..., reverse = FALSE, na.value = "grey85",
                               guide = "colourbar") {
  pal <- andal_pal(reverse)
  ggplot2::scale_fill_gradientn(colours = pal(256), na.value = na.value,
                                guide = guide, ...)
}

scale_color_andal_c <- function(..., reverse = FALSE, na.value = "grey85",
                                guide = "colourbar") {
  pal <- andal_pal(reverse)
  ggplot2::scale_colour_gradientn(colours = pal(256), na.value = na.value,
                                  guide = guide, ...)
}

# Binned / escalonada (tipo scale_fill_viridis_b)
scale_fill_andal_b <- function(..., reverse = FALSE, na.value = "grey85",
                               n.breaks = 6) {
  cols <- c(
    "#052709", "#5FA9B3", "#A17C6C",
    "#D2A899", "#906F00", "#527700"
  )
  if (reverse) cols <- rev(cols)
  ggplot2::scale_fill_stepsn(colors = cols, n.breaks = n.breaks,
                             na.value = na.value, ...)
}

scale_color_andal_b <- function(..., reverse = FALSE, na.value = "grey85",
                                n.breaks = 6) {
  cols <- c(
    "#052709", "#5FA9B3", "#A17C6C",
    "#D2A899", "#906F00", "#527700"
  )
  if (reverse) cols <- rev(cols)
  ggplot2::scale_colour_stepsn(colors = cols, n.breaks = n.breaks,
                               na.value = na.value, ...)
}

# Discreta (por si necesitas categorías)
scale_fill_andal_d <- function(..., reverse = FALSE) {
  cols <- c(
    "#052709", "#5FA9B3", "#A17C6C",
    "#D2A899", "#906F00", "#527700"
  )
  if (reverse) cols <- rev(cols)
  ggplot2::scale_fill_manual(values = cols, ...)
}

scale_color_andal_d <- function(..., reverse = FALSE) {
  cols <- c(
    "#052709", "#5FA9B3", "#A17C6C",
    "#D2A899", "#906F00", "#527700"
  )
  if (reverse) cols <- rev(cols)
  ggplot2::scale_colour_manual(values = cols, ...)
}


Estructurar los datos

Code
# Function to clean and prepare region names for matching
# Load the fuzzy matching library
library(stringdist)

Adjuntando el paquete: ‘stringdist’

The following object is masked from ‘package:terra’:

extract

The following object is masked from ‘package:tidyr’:

extract
Code
ref_regions <- tibble::tribble(
  ~roman_code, ~canon_name, ~aliases,
  "XV",  "arica y parinacota",                    c("arica y parinacota","arica","parinacota"),
  "I",   "tarapaca",                               c("tarapaca"),
  "II",  "antofagasta",                            c("antofagasta"),
  "III", "atacama",                                c("atacama"),
  "IV",  "coquimbo",                               c("coquimbo"),
  "V",   "valparaiso",                             c("valparaíso","valparaiso"),
  "RM",  "metropolitana",                          c("rm","metropolitana","metropolitana de santiago","santiago"),
  "VI",  "ohiggins",                               c("o'higgins","ohiggins","libertador general bernardo ohiggins","bernardo ohiggins","libertador ohiggins"),
  "VII", "maule",                                  c("maule"),
  "VIII","biobio",                                 c("biobio","bio bio","bio-bio","bío-bío"),
  "IX",  "araucania",                              c("la araucania","araucania","araucanía"),
  "X",   "los lagos",                              c("los lagos"),
  "XI",  "aysen",                                  c("aysen","aysén"),
  "XII", "magallanes y la antartica chilena",      c("magallanes","antartica","antártica","magallanes y antartica","magallanes y la antartica chilena"),
  "XIV", "los rios",                               c("los rios","los ríos")
)

# --- Limpieza homogénea (preserva 'los/las') ---
clean_string <- function(s){
  s <- tolower(s)
  s <- iconv(s, from = "UTF-8", to = "ASCII//TRANSLIT")  # quita tildes
  s <- gsub("[[:punct:]]", " ", s)                        # deja espacios
  s <- gsub("[[:digit:]]", "", s)
  # quitamos conectores pero NO 'los/las' (desambiguación)
  s <- gsub("\\b(de|del|y|la|el|region|región)\\b", " ", s)
  s <- gsub("\\s+", " ", s)
  trimws(s)
}

# expandimos alias y limpiaremos todos con la misma función
ref_alias <- ref_regions |>
  tidyr::unnest_longer(aliases) |>
  dplyr::mutate(clean_alias = clean_string(aliases),
                clean_canon = clean_string(canon_name)) |>
  dplyr::distinct(roman_code, canon_name, clean_alias, clean_canon)

# --- Reglas exactas/regex muy específicas (antes de fuzzy) ---
rule_match_code <- function(x_clean){
  # pares conflictivos: exigir frase completa
  if(grepl("\\blos\\s+lagos\\b", x_clean)) return("X")
  if(grepl("\\blos\\s+rios\\b",  x_clean)) return("XIV")

  # otros casos robustos
  if(grepl("\\bbiobio\\b|\\bbio\\s*bio\\b", x_clean)) return("VIII")
  if(grepl("\\bo\\s*higgins\\b|\\blibertador\\b|\\bbernardo\\b", x_clean)) return("VI")
  if(grepl("\\bmetropolitana\\b|\\brm\\b|\\bsantiago\\b", x_clean)) return("RM")
  if(grepl("\\bmagallanes\\b|\\bantartic", x_clean)) return("XII")
  if(grepl("\\baraucan", x_clean)) return("IX")
  if(grepl("\\bvalparais", x_clean)) return("V")
  if(grepl("\\btarapac", x_clean)) return("I")
  if(grepl("\\barica\\b|\\bparinacota\\b", x_clean)) return("XV")
  if(grepl("\\bantofagasta\\b", x_clean)) return("II")
  if(grepl("\\batacama\\b", x_clean)) return("III")
  if(grepl("\\bcoquimbo\\b", x_clean)) return("IV")
  if(grepl("\\bmaule\\b", x_clean)) return("VII")
  if(grepl("\\baysen\\b", x_clean)) return("XI")
  return(NA_character_)
}

# --- Matcher seguro ---
match_regions_safe <- function(messy_names, ref_alias,
                               jw_thresh = 0.12, tie_margin = 0.02){

  sapply(messy_names, function(x){
    x_clean <- clean_string(x)

    # 1) Reglas
    code <- rule_match_code(x_clean)
    if(!is.na(code)) return(code)

    # 2) Igualdad exacta con alias limpios o canónicos
    eq_hit <- ref_alias |>
      dplyr::filter(x_clean == clean_alias | x_clean == clean_canon)
    if(nrow(eq_hit) == 1) return(eq_hit$roman_code[[1]])
    if(nrow(eq_hit) > 1)  return(NA_character_)  # ambigüedad improbable: revisar

    # 3) Fuzzy con umbral y desempate
    cand <- unique(c(ref_alias$clean_alias, ref_alias$clean_canon))
    d <- stringdist::stringdist(x_clean, cand, method = "jw")
    i_min <- which.min(d); d_min <- d[i_min]

    # si la mejor distancia es mala, no arriesgar
    if(is.infinite(d_min) || is.na(d_min) || d_min > jw_thresh) return(NA_character_)

    # empate cercano: si la diferencia entre 1º y 2º es pequeña, exigir frase completa o marcar NA
    ord <- order(d); if(length(ord) >= 2){
      if((d[ord[2]] - d[ord[1]]) < tie_margin) return(NA_character_)
    }

    cand_str <- cand[i_min]
    hit <- ref_alias |>
      dplyr::filter(clean_alias == cand_str | clean_canon == cand_str) |>
      dplyr::slice(1)
    hit$roman_code[[1]]
  })
}

# This is your sample data from the query


matched <- match_regions_safe(CONS_C3_mod_c_a2$region_del_centro, ref_alias)

# Add the new codes back to your original data frame
CONS_C3_mod_c_a2$CC1 <- matched

table(CONS_C3_mod_c_a2$region_del_centro, CONS_C3_mod_c_a2$CC1) |> 
  knitr::kable("markdown", caption="Tabla de contingencia entre nombres originales y códigos asignados")


# Now you can proceed with your lexpand call, including 'reg_res' in the aggre list
SISTRAT_c3_reg <- popEpi::lexpand(
  CONS_C3_mod_c_a2,
  status = status,
  birth = birth_date_rec_joined,
  exit = def_date_na,
  entry = fecha_ingreso_a_tratamiento,
  # 'reg_res' has been REMOVED from the breaks list
  breaks = list(
    per = seq(2011, 2021, by = 1), #En popEpi::lexpand, los breaks son cortes y generan intervalos [a, b) (cerrado a la izquierda, abierto a la derecha). Cambio en 13-octubre 2025 2025-10-13
    age = c(18, 30, 45, 60, Inf)
  ),
  # 'reg_res' has been ADDED to the aggre list
  aggre = list(
    agegroup = age,
    sex = sex,
    CC1 = CC1
  )
)

# --- 1) Prepare SISTRAT cohort by CC_1 × sex × age ---
# from0to1 is deaths; pyrs are person-years
SISTRAT_agg <- SISTRAT_c3_reg %>%
  mutate(
    sex_rec = recode(tolower(sex), male = "hombre", female = "mujer"),
    adm_age_cat = case_when(
      agegroup == 18 ~ "18-29",
      agegroup == 30 ~ "30-44",
      agegroup == 45 ~ "45-59",
      agegroup == 60 ~ "60-78",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(adm_age_cat), !is.na(CC1)) %>%
  group_by(CC1, sex_rec, adm_age_cat) %>%
  summarise(
    Y  = sum(from0to1, na.rm = TRUE),
    PY = sum(pyrs,      na.rm = TRUE),
    .groups = "drop"
  )

# Cohort totals by CC_1 (for modeling)
SISTRAT_by_CC1 <- SISTRAT_agg %>%
  group_by(CC1) %>%
  summarise(Y = sum(Y), PY = sum(PY), .groups = "drop")
Tabla de contingencia entre nombres originales y códigos asignados
I II III IV IX RM V VI VII VIII X XII XIV XV
de antofagasta 0 119 0 0 0 0 0 0 0 0 0 0 0 0
de arica y parinacota 0 0 0 0 0 0 0 0 0 0 0 0 0 102
de atacama 0 0 22 0 0 0 0 0 0 0 0 0 0 0
de coquimbo 0 0 0 68 0 0 0 0 0 0 0 0 0 0
de la araucania 0 0 0 0 42 0 0 0 0 0 0 0 0 0
de los lagos 0 0 0 0 0 0 0 0 0 0 81 0 0 0
de los rios 0 0 0 0 0 0 0 0 0 0 0 0 56 0
de magallanes y la antartica chilena 0 0 0 0 0 0 0 0 0 0 0 78 0 0
de tarapaca 29 0 0 0 0 0 0 0 0 0 0 0 0 0
de valparaiso 0 0 0 0 0 0 196 0 0 0 0 0 0 0
del bio-bio 0 0 0 0 0 0 0 0 0 154 0 0 0 0
del libertador general bernardo ohiggins 0 0 0 0 0 0 0 82 0 0 0 0 0 0
del maule 0 0 0 0 0 0 0 0 41 0 0 0 0 0
metropolitana 0 0 0 0 0 321 0 0 0 0 0 0 0 0

Definir la población de referencia y su mortalidad

Code
# Load the reference mortality rates
proy_ine_reg_group_25_oct_reg_al2015_dput<- 
structure(list(reg_res = c("01", "01", "01", "01", "01", "01", 
"01", "01", "02", "02", "02", "02", "02", "02", "02", "02", "03", 
"03", "03", "03", "03", "03", "03", "03", "04", "04", "04", "04", 
"04", "04", "04", "04", "05", "05", "05", "05", "05", "05", "05", 
"05", "06", "06", "06", "06", "06", "06", "06", "06", "07", "07", 
"07", "07", "07", "07", "07", "07", "08", "08", "08", "08", "08", 
"08", "08", "08", "09", "09", "09", "09", "09", "09", "09", "09", 
"10", "10", "10", "10", "10", "10", "10", "10", "11", "11", "11", 
"11", "11", "11", "11", "11", "12", "12", "12", "12", "12", "12", 
"12", "12", "13", "13", "13", "13", "13", "13", "13", "13", "14", 
"14", "14", "14", "14", "14", "14", "14", "15", "15", "15", "15", 
"15", "15", "15", "15", "16", "16", "16", "16", "16", "16", "16", 
"16"), year = c(2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015), 
    sex = c("Male", "Male", "Male", "Male", "Female", "Female", 
    "Female", "Female", "Male", "Male", "Male", "Male", "Female", 
    "Female", "Female", "Female", "Male", "Male", "Male", "Male", 
    "Female", "Female", "Female", "Female", "Male", "Male", "Male", 
    "Male", "Female", "Female", "Female", "Female", "Male", "Male", 
    "Male", "Male", "Female", "Female", "Female", "Female", "Male", 
    "Male", "Male", "Male", "Female", "Female", "Female", "Female", 
    "Male", "Male", "Male", "Male", "Female", "Female", "Female", 
    "Female", "Male", "Male", "Male", "Male", "Female", "Female", 
    "Female", "Female", "Male", "Male", "Male", "Male", "Female", 
    "Female", "Female", "Female", "Male", "Male", "Male", "Male", 
    "Female", "Female", "Female", "Female", "Male", "Male", "Male", 
    "Male", "Female", "Female", "Female", "Female", "Male", "Male", 
    "Male", "Male", "Female", "Female", "Female", "Female", "Male", 
    "Male", "Male", "Male", "Female", "Female", "Female", "Female", 
    "Male", "Male", "Male", "Male", "Female", "Female", "Female", 
    "Female", "Male", "Male", "Male", "Male", "Female", "Female", 
    "Female", "Female", "Male", "Male", "Male", "Male", "Female", 
    "Female", "Female", "Female"), agegroup = c(18, 30, 45, 60, 
    18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 
    60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 
    45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 
    30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 
    18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 
    60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 
    45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 
    30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 18, 30, 45, 60, 
    18, 30, 45, 60), poblacion = c(34002L, 38653L, 28302L, 15091L, 
    32200L, 36648L, 27375L, 16273L, 64007L, 72673L, 53668L, 27557L, 
    60944L, 67667L, 51267L, 30511L, 28880L, 32559L, 28056L, 16330L, 
    27350L, 30400L, 27570L, 17325L, 71827L, 78583L, 69004L, 45243L, 
    71480L, 81367L, 72797L, 51981L, 180134L, 186088L, 167910L, 
    118195L, 172861L, 186938L, 182401L, 142333L, 79763L, 101327L, 
    93062L, 58682L, 76880L, 99496L, 92575L, 63692L, 93197L, 110899L, 
    103897L, 68655L, 96053L, 113384L, 106240L, 75650L, 155816L, 
    164863L, 152554L, 94068L, 155474L, 171681L, 161491L, 111956L, 
    88721L, 99564L, 90612L, 59966L, 91078L, 103766L, 92400L, 
    68850L, 77395L, 96278L, 82763L, 48381L, 75995L, 95497L, 80050L, 
    54458L, 9138L, 12044L, 10142L, 5664L, 8611L, 11949L, 9318L, 
    5554L, 15950L, 19178L, 16790L, 10861L, 14332L, 18242L, 16390L, 
    10938L, 760210L, 819336L, 644988L, 385795L, 727290L, 802028L, 
    699815L, 487007L, 37117L, 38145L, 39261L, 24363L, 37110L, 
    39544L, 38848L, 28148L, 23852L, 24685L, 19338L, 12961L, 21661L, 
    24452L, 20876L, 15206L, 41936L, 47039L, 49239L, 33580L, 43939L, 
    50966L, 51433L, 37434L)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -128L))
#mort_2015_reg
mort_2015_reg_dput<- 
structure(list(reg_res = c("01", "01", "01", "01", "01", "01", 
"01", "01", "01", "01", "02", "02", "02", "02", "02", "02", "02", 
"02", "02", "02", "03", "03", "03", "03", "03", "03", "03", "03", 
"03", "03", "04", "04", "04", "04", "04", "04", "04", "04", "04", 
"04", "05", "05", "05", "05", "05", "05", "05", "05", "05", "05", 
"06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "07", 
"07", "07", "07", "07", "07", "07", "07", "07", "07", "07", "08", 
"08", "08", "08", "08", "08", "08", "08", "08", "08", "08", "09", 
"09", "09", "09", "09", "09", "09", "09", "09", "09", "09", "10", 
"10", "10", "10", "10", "10", "10", "10", "10", "10", "10", "11", 
"11", "11", "11", "11", "11", "11", "11", "11", "11", "12", "12", 
"12", "12", "12", "12", "12", "12", "12", "12", "13", "13", "13", 
"13", "13", "13", "13", "13", "13", "13", "13", "14", "14", "14", 
"14", "14", "14", "14", "14", "14", "14", "14", "15", "15", "15", 
"15", "15", "15", "15", "15", "15", "15", "15"), sexo = c(1, 
1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 
1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 
1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 
1, 2, 2, 2, 2, 2, 9, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 9, 1, 1, 1, 
1, 1, 2, 2, 2, 2, 2, 9, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 9, 1, 1, 
1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1, 
1, 1, 2, 2, 2, 2, 2, 9, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 9, 1, 1, 
1, 1, 1, 2, 2, 2, 2, 2, 9), adm_age_cat = c("18-29", "30-44", 
"45-59", "60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, 
"18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", "45-59", 
"60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", 
"30-44", "45-59", "60-78", NA, "18-29", "30-44", "45-59", "60-78", 
NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", 
"45-59", "60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, 
"18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", "45-59", 
"60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", 
"30-44", "45-59", "60-78", NA, NA, "18-29", "30-44", "45-59", 
"60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, NA, "18-29", 
"30-44", "45-59", "60-78", NA, "18-29", "30-44", "45-59", "60-78", 
NA, NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", 
"45-59", "60-78", NA, NA, "18-29", "30-44", "45-59", "60-78", 
NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", 
"45-59", "60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, 
"18-29", "30-44", "45-59", "60-78", NA, "18-29", "30-44", "45-59", 
"60-78", NA, NA, "18-29", "30-44", "45-59", "60-78", NA, "18-29", 
"30-44", "45-59", "60-78", NA, NA, "18-29", "30-44", "45-59", 
"60-78", NA, "18-29", "30-44", "45-59", "60-78", NA, NA), n = c(40L, 
73L, 134L, 316L, 228L, 11L, 31L, 57L, 196L, 325L, 71L, 99L, 338L, 
729L, 483L, 28L, 47L, 144L, 477L, 617L, 23L, 69L, 137L, 321L, 
326L, 17L, 20L, 73L, 202L, 366L, 83L, 98L, 355L, 895L, 878L, 
21L, 61L, 203L, 561L, 1001L, 179L, 271L, 811L, 2531L, 2516L, 
45L, 148L, 536L, 1824L, 3413L, 97L, 172L, 449L, 1167L, 1072L, 
33L, 95L, 264L, 795L, 1287L, 120L, 184L, 494L, 1475L, 1326L, 
43L, 92L, 292L, 978L, 1485L, 2L, 188L, 384L, 1127L, 2850L, 2462L, 
71L, 153L, 667L, 1953L, 3057L, 1L, 98L, 206L, 506L, 1359L, 1331L, 
33L, 93L, 303L, 938L, 1631L, 4L, 101L, 202L, 507L, 1026L, 952L, 
31L, 82L, 229L, 748L, 1228L, 2L, 15L, 22L, 64L, 119L, 99L, 1L, 
7L, 21L, 61L, 89L, 18L, 28L, 89L, 247L, 176L, 6L, 16L, 45L, 164L, 
244L, 700L, 1169L, 3089L, 7822L, 6991L, 238L, 565L, 1873L, 5714L, 
10914L, 8L, 42L, 96L, 215L, 583L, 549L, 14L, 35L, 106L, 383L, 
566L, 1L, 33L, 42L, 99L, 290L, 233L, 8L, 23L, 52L, 182L, 282L, 
1L), sex_rec = c("hombre", "hombre", "hombre", "hombre", "hombre", 
"mujer", "mujer", "mujer", "mujer", "mujer", "hombre", "hombre", 
"hombre", "hombre", "hombre", "mujer", "mujer", "mujer", "mujer", 
"mujer", "hombre", "hombre", "hombre", "hombre", "hombre", "mujer", 
"mujer", "mujer", "mujer", "mujer", "hombre", "hombre", "hombre", 
"hombre", "hombre", "mujer", "mujer", "mujer", "mujer", "mujer", 
"hombre", "hombre", "hombre", "hombre", "hombre", "mujer", "mujer", 
"mujer", "mujer", "mujer", "hombre", "hombre", "hombre", "hombre", 
"hombre", "mujer", "mujer", "mujer", "mujer", "mujer", "hombre", 
"hombre", "hombre", "hombre", "hombre", "mujer", "mujer", "mujer", 
"mujer", "mujer", "mujer", "hombre", "hombre", "hombre", "hombre", 
"hombre", "mujer", "mujer", "mujer", "mujer", "mujer", "mujer", 
"hombre", "hombre", "hombre", "hombre", "hombre", "mujer", "mujer", 
"mujer", "mujer", "mujer", "mujer", "hombre", "hombre", "hombre", 
"hombre", "hombre", "mujer", "mujer", "mujer", "mujer", "mujer", 
"mujer", "hombre", "hombre", "hombre", "hombre", "hombre", "mujer", 
"mujer", "mujer", "mujer", "mujer", "hombre", "hombre", "hombre", 
"hombre", "hombre", "mujer", "mujer", "mujer", "mujer", "mujer", 
"hombre", "hombre", "hombre", "hombre", "hombre", "mujer", "mujer", 
"mujer", "mujer", "mujer", "mujer", "hombre", "hombre", "hombre", 
"hombre", "hombre", "mujer", "mujer", "mujer", "mujer", "mujer", 
"mujer", "hombre", "hombre", "hombre", "hombre", "hombre", "mujer", 
"mujer", "mujer", "mujer", "mujer", "mujer")), row.names = c(NA, 
-157L), class = c("tbl_df", "tbl", "data.frame")) 

mort_2015_reg_dput<- 
dplyr::filter(data.frame(mort_2015_reg_dput),!is.na(adm_age_cat))

# Mapeo de número romano (CC_1) -> código CUT (reg_res)
roman_to_cut <- tibble::tibble(
  CC_1 = c("I","II","III","IV","V","VI","VII","VIII","IX","X","XI","XII","RM","XIV","XV","XVI"),
  reg_res = c("01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16")
)


# --- 2) Build background reference rates by CC_1 × sex × age (2015) ---
# Population
pop <- proy_ine_reg_group_25_oct_reg_al2015_dput|>
  filter(year == 2015, sex %in% c("Male","Female"))|>
  dplyr::mutate(
    sex_rec = recode(sex, Male = "hombre", Female = "mujer"),
    adm_age_cat = case_when(
      agegroup == 18 ~ "18-29",
      agegroup == 30 ~ "30-44",
      agegroup == 45 ~ "45-59",
      agegroup == 60 ~ "60-78",
      TRUE ~ NA_character_
    ),
    reg_res = stringr::str_pad(reg_res, width = 2, pad = "0")
  )|>
  dplyr::filter(!is.na(adm_age_cat))|>
  dplyr::select(reg_res, sex_rec, adm_age_cat, poblacion)|>
  dplyr::left_join(roman_to_cut, by = "reg_res")|> 
  dplyr::rename("CC1" = "CC_1")

# Deaths
deaths <- mort_2015_reg_dput|>
  dplyr::filter(!is.na(adm_age_cat), sexo %in% c(1,2))|>
  dplyr::mutate(
    sex_rec = if_else(sexo == 1, "hombre", "mujer"),
    reg_res = stringr::str_pad(reg_res, width = 2, pad = "0"))|>
  dplyr::select(reg_res, sex_rec, adm_age_cat, n)|>
  dplyr::left_join(roman_to_cut, by = "reg_res")|> 
  dplyr::rename("CC1" = "CC_1")

# Make complete shells to avoid missing strata
shell_cc1 <- expand.grid(
  CC1       = unique(roman_to_cut$CC_1),
  sex_rec    = c("hombre","mujer"),
  adm_age_cat= c("18-29","30-44","45-59","60-78"),
  stringsAsFactors = FALSE
) |> tibble::as_tibble()

pop_ref <- shell_cc1|>
  left_join(pop,    by = c("CC1","sex_rec","adm_age_cat"))|>
  mutate(poblacion = replace_na(poblacion, 0L))

deaths_ref <- shell_cc1|>
  left_join(deaths, by = c("CC1","sex_rec","adm_age_cat"))|>
  mutate(n = replace_na(n, 0L))

# Reference rates r_{CC_1,a,s} = D / P
P_df <- pop_ref  |> group_by(CC1, sex_rec, adm_age_cat)|> summarise(P = sum(poblacion), .groups = "drop")
D_df <- deaths_ref%>% group_by(CC1, sex_rec, adm_age_cat)|> summarise(D = sum(n),         .groups = "drop")
rates_ref <- full_join(D_df, P_df, by = c("CC1","sex_rec","adm_age_cat"))|>
  mutate(
    D = coalesce(D, 0L),
    P = coalesce(P, 0L),
    rate = if_else(P > 0, D / P, 0)
  )

# --- 3) Expected deaths for the cohort by CC_1: E = Σ PY × rate_ref ---
E_by_CC1 <- SISTRAT_agg|>
  left_join(rates_ref, by = c("CC1","sex_rec","adm_age_cat"))|>
  mutate(rate = replace_na(rate, 0))|>
  group_by(CC1)|>
  summarise(
    E  = sum(PY * rate),
    PY = sum(PY),
    Y  = sum(Y),
    .groups = "drop"
  )|>
  mutate(E = if_else(E == 0, 1e-8, E))   # guard against exact zeros

Formatear el mapa y matriz de adyacencia

Armonizamos el mapa con los datos y construimos la matriz de adyacencia para el modelo BYM, fijando a priori conservadores: asumimos una probabilidad baja (≈1%) de que la desviación estándar del efecto aleatorio supere 1 (con esto, asumimos poca variabilidad interregional) no y predefinimos un equilibrio entre el componente espacial estructurado (por vecindad) y el no estructurado (ruido, o particularidad de la región).

Code
# --- 4) Geometry (join by CC_1) and adjacency for BYM ---
chile_sf <- geodata::gadm(country = "CHL", level = 1, path = tempdir()) |>
  sf::st_as_sf() |>
  dplyr::mutate(CC_1 = as.character(CC_1))|>
  # keep only regions we can map (GADM may or may not have XVI; most recent does)
  (\(df) {
    if (filter(df, !CC_1 %in% roman_to_cut$CC_1) |> nrow()>0) {
      message("Note: data does not include X10=Región de Los Lagos; XI11=Región de Aysén del General Carlos Ibáñez del Campo; XVI16=Región de Ñuble")
      print(filter(df, !CC_1 %in% roman_to_cut$CC_1) |>pull(CC1))
    }
    df
  })()|>
  dplyr::filter(CC_1 %in% roman_to_cut$CC_1)|>
  dplyr::arrange(CC_1)|>
  dplyr::mutate(region_id = dplyr::row_number())

# Align modeling data to shapes (right_join keeps all regions in the map)
bym_df <- E_by_CC1|>
  dplyr::right_join(sf::st_drop_geometry(chile_sf)|> dplyr::select(CC_1, region_id), by = c("CC1"="CC_1"))|>
  dplyr::arrange(region_id)|>
  dplyr::mutate(
    Y  = replace_na(Y, 0),
    PY = replace_na(PY, 0),
    E  = replace_na(E,  1e-8)  # tiny to allow model to run if truly no expected deaths
  )

# Adjacency
nb <- spdep::poly2nb(chile_sf, queen = TRUE)
adj_file <- file.path(tempdir(), "chile_cc1_adj.graph")
spdep::nb2INLA(file = adj_file, nb)

# Priors
hyper <- list(
  prec = list(prior = "pc.prec", param = c(1, 0.01)),
  phi  = list(prior = "pc", param = c(0.5, 2/3))
)

Modelo Bayesiano BYM2 por tasas suavizadas, riesgo relativo y muertes en exceso

Obtenemos las tasas por 100 años-persona, riesgo relativo y muertes en exceso suavizadas por vecindad espacial, además de la probabilidad de que el riesgo relativo sea mayor que 1.

Code
# --- 5A) BYM for smoothed death rates (cohort) ---
# Y ~ Poisson(PY × rate_i); fitted values are on scale of rate (per person-year)
fit_rate <- INLA::inla(
  Y ~ 1 + f(region_id, model = "bym2", graph = adj_file, scale.model = TRUE, hyper = hyper),
  data = bym_df,
  family = "poisson",
  E = PY,  # exposure = person-years
  control.predictor = list(compute = TRUE),
  control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE)
)

# Smoothed rates (per person-year)
rate_mean <- fit_rate$summary.fitted.values$mean      # rate
rate_lcl  <- fit_rate$summary.fitted.values$`0.025quant`
rate_ucl  <- fit_rate$summary.fitted.values$`0.975quant`

# If you prefer per 1,000 PY:
per   <- 100
rate_mean_100 <- rate_mean * per
rate_lcl_100  <- rate_lcl  * per
rate_ucl_100  <- rate_ucl  * per

# --- 5B) BYM for excess risk and excess deaths ---
# Y ~ Poisson(E × RR_i); fitted values are RR
fit_rr <- INLA::inla(
  Y ~ 1 + f(region_id, model = "bym2", graph = adj_file, scale.model = TRUE, hyper = hyper),
  data = bym_df,
  family = "poisson",
  E = E,   # expected deaths from reference
  control.predictor = list(compute = TRUE),
  control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE)
)

RR_mean <- fit_rr$summary.fitted.values$mean
RR_lcl  <- fit_rr$summary.fitted.values$`0.025quant`
RR_ucl  <- fit_rr$summary.fitted.values$`0.975quant`

# Smoothed excess deaths = (RR - 1) * E
Excess_mean <- (RR_mean - 1) * bym_df$E
Excess_lcl  <- (RR_lcl  - 1) * bym_df$E
Excess_ucl  <- (RR_ucl  - 1) * bym_df$E

# Exceedance probability: P(RR > 1) = P(eta > 0)
if (!is.null(fit_rr$marginals.linear.predictor) &&
    length(fit_rr$marginals.linear.predictor) == nrow(bym_df)) {
  p_exceed_1 <- sapply(fit_rr$marginals.linear.predictor, function(m) 1 - INLA::inla.pmarginal(0, m))
} else {
  # Monte Carlo backup
  smp <- INLA::inla.posterior.sample(4000, fit_rr)
  pred_names <- grep("^Predictor", rownames(smp[[1]]$latent), value = TRUE)
  lp_mat <- sapply(smp, function(s) as.numeric(s$latent[pred_names]))
  p_exceed_1 <- rowMeans(lp_mat > 0, na.rm = TRUE)
}


ps_pois <- inla.posterior.sample(2000, fit_rr, selection = list(Predictor = 1:nrow(fit_rr$summary.linear.predictor)))
eta_mat_pois <- sapply(ps_pois, function(s) s$latent[grep("^Predictor", rownames(s$latent)), , drop = TRUE])
p_exceed_1_mc <- rowMeans(eta_mat_pois > 0)

# --- 6) Bind results and map ---
res_cc1 <- bym_df|>
  mutate(
    rate_mean      = rate_mean,
    rate_lcl       = rate_lcl,
    rate_ucl       = rate_ucl,
    rate_mean_100 = rate_mean_100,
    rate_lcl_100  = rate_lcl_100,
    rate_ucl_100  = rate_ucl_100,
    RR_mean        = RR_mean,
    RR_lcl         = RR_lcl,
    RR_ucl         = RR_ucl,
    Excess_mean    = Excess_mean,
    Excess_lcl     = Excess_lcl,
    Excess_ucl     = Excess_ucl,
    p_exceed_1     = p_exceed_1_mc
  )

result_sf <- chile_sf|>
  left_join(res_cc1, by = c("CC_1"="CC1"))

orden_geografico_romano <- c(
  "XV",   # Arica y Parinacota
  "I",    # Tarapacá
  "II",   # Antofagasta
  "III",  # Atacama
  "IV",   # Coquimbo
  "V",    # Valparaíso
  "RM",   # Metropolitana
  "VI",   # O'Higgins
  "VII",  # Maule
  "XVI",  # Ñuble
  "VIII", # Biobío
  "IX",   # La Araucanía
  "XIV",  # Los Ríos
  "X",    # Los Lagos
  "XI",   # Aysén (lo incluyo por si aparece en tus datos)
  "XII"   # Magallanes
)


res_cc1 |> 
  dplyr::select(-dplyr::any_of(c("rate_mean","rate_lcl","rate_ucl", "p_exceed_1")))|>
   dplyr::mutate(CC1 = factor(CC1, levels = orden_geografico_romano),  )|>
  dplyr::arrange(CC1)|>
  knitr::kable("markdown", digits=1, caption= "Resultados del modelo BYM2 por Región")
Resultados del modelo BYM2 por Región
CC1 E PY Y region_id rate_mean_100 rate_lcl_100 rate_ucl_100 RR_mean RR_lcl RR_ucl Excess_mean Excess_lcl Excess_ucl
XV 1.8 427.2 7 15 2.1 1.2 3.4 6.0 3.4 8.5 8.9 4.4 13.4
I 0.8 186.5 3 1 2.4 1.2 4.1 6.4 3.7 9.3 4.3 2.1 6.6
II 1.6 413.9 12 2 2.9 1.8 4.4 7.2 4.9 10.2 10.0 6.3 14.9
III 0.4 96.1 9 3 4.7 2.5 8.9 9.1 5.8 16.7 3.3 2.0 6.4
IV 0.8 200.5 8 4 3.4 2.1 5.6 7.7 5.2 11.8 5.5 3.4 8.8
V 3.2 765.6 17 7 2.5 1.7 3.5 6.4 4.4 8.5 17.3 10.8 23.8
RM 3.9 1251.5 32 6 2.6 1.9 3.5 7.6 5.8 10.1 25.6 18.5 35.2
VI 1.3 306.3 5 8 2.4 1.4 3.8 6.3 3.8 8.7 7.1 3.7 10.4
VII 1.0 220.3 6 9 3.0 1.7 4.9 6.9 4.4 9.9 5.9 3.4 8.9
XVI 0.0 0.0 0 16 3.7 1.7 7.4 7.3 4.2 11.8 0.0 0.0 0.0
VIII 3.3 551.2 19 10 3.7 2.5 5.3 6.7 4.6 8.8 18.8 12.0 25.8
IX 0.7 76.0 7 5 5.9 3.3 10.5 7.9 5.3 12.2 5.0 3.1 8.0
XIV 1.3 195.2 11 14 5.6 3.4 8.5 7.9 5.4 11.5 9.1 5.9 13.9
X 1.6 219.6 18 11 7.0 4.5 10.5 9.0 6.3 13.5 12.4 8.2 19.5
XI 0.0 0.0 0 12 6.3 2.6 12.6 8.1 4.6 13.9 0.0 0.0 0.0
XII 2.1 276.6 18 13 6.3 4.1 9.4 8.1 5.8 11.6 15.0 10.0 22.4

Visualización de resultados

Code
# Print plots
cowplot::plot_grid(
    p_tratados +
        coord_sf(expand = FALSE) +
        theme(
            plot.margin = margin(0, 0, 0, 0),
            legend.spacing = unit(0.1, "lines"),
            legend.key.size = unit(0.6, "lines")
        ), p_deaths +
        coord_sf(expand = FALSE) +
        theme(
            plot.margin = margin(0, 0, 0, 0),
            legend.spacing = unit(0.1, "lines"),
            legend.key.size = unit(0.6, "lines")
        ),
    ggplot(result_sf) +
        geom_sf(aes(fill = RR_mean), color = "grey60", size = 0.2) +
        scale_fill_gradientn(colors = rev(c(g052709, dark_teal, c5fa9b3, d_edee4, blanco)), na.value = "#DEDEE4", name = "RME Suavizado") +
        #scale_fill_viridis_c(option = "C", name = "Smoothed RR") +
        labs(title = NULL,#"BYM2-smoothed excess risk (SISTRAT vs reference)",
             subtitle = NULL)+#"RR = observed / expected (age×sex region-specific rates)") +
        theme_void(base_size = 14) +
        #xlim(c(80,66.4))+
        xlim(c(-75, -66.4))+
        theme(legend.position = "left",
              legend.title = element_text(size=12.9),
              legend.text = element_text(size=11.9))+
        ggrepel::geom_label_repel(
            data =  dplyr::mutate(result_sf, 
                                  # 1. Calculates a point guaranteed to be within the polygon boundaries.
                                  point_on_surface = sf::st_point_on_surface(geometry),
                                  
                                  # 2. Extracts the X (Longitude) and Y (Latitude) coordinates from that point.
                                  lon = sf::st_coordinates(point_on_surface)[, "X"],
                                  lat = sf::st_coordinates(point_on_surface)[, "Y"]), # Usamos los datos con las coordenadas calculadas
            aes(x = lon, y = lat, label = CC_1), # Coordenadas y texto de la etiqueta
            size = 4.6,                              # Tamaño del texto
            min.segment.length = 0.5,                # Longitud mínima de la línea de conexión
            box.padding = 0.3,                       # Espacio alrededor del texto dentro de la caja
            point.padding = 0.5,                     # Distancia de la caja a la geometría
            # Fondo de la etiqueta
            fill = alpha("white", 0.7),              # Color de fondo semi-transparente
            max.overlaps = Inf,    # Increase repulsion effort to avoid all overlaps
            seed = 2125,             # For reproducible results
            color = "black"                          # Color del texto
        ),
    nrow = 1,                # o nrow = 1 según tu disposición
    align = "v",             # o "h" para horizontal
    rel_heights = c(1, 1),# ajusta proporciones si necesario
    axis = "tblr",           # alinea ejes si aplica
    labels = NULL,           # evita espacio extra por etiquetas
    label_size = 0,          # si usas labels, reduce tamaño
    label_x = 0, label_y = 1,# controla posición de etiquetas
    vjust = .5, hjust = 0,    # controla alineación vertical/horizontal
    scale = c(0.99, .99, 1.03) # reduce tamaño relativo de cada plot
)

Coordinate system already present. ℹ Adding new coordinate system, which will replace the existing one. Coordinate system already present. ℹ Adding new coordinate system, which will replace the existing one.

Warning: There was 1 warning in stopifnot(). ℹ In argument: point_on_surface = sf::st_point_on_surface(geometry). Caused by warning in st_point_on_surface.sfc(): ! st_point_on_surface may not give correct results for longitude/latitude data

Code
if(grepl("ubuntu",pak::system_r_platform())){
  ggsave(paste0("map_con_RR.png"), dpi = 600, scale=1,   width = 12, height = 7.38)
} else{
  ggsave(paste0("map_con_RR.png"), dpi = 600, scale=1,   width = 12, height = 7.38)
}

Resumen

Code
# 1) Orden final (Totales arriba). Elimina la línea con `sections_order <- rev(sections_order)`
sections_order <- c(
  "Totales",
  "Por causas (aj. por sexo y edad)",
  "Edad (aj. por sexo)",
  "Sexo (aj. por grupo de edad)"
)

capitalize_first_letter_robust <- function(x) {
  # 1. Limpiar espacios al inicio y final
  x_clean <- trimws(x)

  # 2. Capitalizar la primera letra
  first_letter <- toupper(substring(x_clean, 1, 1))
  
  # 3. Tomar el resto de la cadena a partir del segundo caracter
  rest_of_string <- substring(x_clean, 2)

  return(paste0(first_letter, rest_of_string))
}


# --- (después de tu asignación manual de `smr_resumen$section <- c(...)`) ---

# 3) Mantener solo esas secciones y fijar el orden como factor
keep_sections <- sections_order
smr_resumen2 <- smr_resumen |>
  dplyr::filter(section %in% keep_sections) |>
  dplyr::mutate(section = factor(section, levels = sections_order))|> 
  dplyr::filter(!res %in% c("Total ajustado por grupos de edad, sexo y año calendario", "Total ajustado por edad, sexo y año calendario"))|> 
  dplyr::mutate(section= dplyr::case_when(grepl("causas",section)~"Por causas (aj. por sexo\ny grupos de edad)",
                                           TRUE~section))|> 
  dplyr::mutate(section= gsub("grupo de","grupos de", section))|> 
  dplyr::mutate(label = gsub("^Total", "", label))|>
  dplyr::mutate(label = sapply(label, capitalize_first_letter_robust))


sections_order <- c(
  "Totales",
  "Por causas (aj. por sexo\ny grupos de edad)",
  "Edad (aj. por sexo)",
  "Sexo (aj. por grupos de edad)"
)

# Create display
display_list <- list()
y <- 0
for (sec in sections_order) {
  header <- tibble(section = sec, label = NA_character_, obs = NA_real_, py = NA_real_, exp = NA_real_,
                   rme = NA_real_, lo = NA_real_, hi = NA_real_, is_header = TRUE, y = y)
  display_list <- append(display_list, list(header))
  y <- y + 1

  items <- smr_resumen2 %>% filter(section == sec)
  if (nrow(items) > 0) {
    items <- items %>% mutate(is_header = FALSE, y = seq(y, by = 1, length.out = n()))
    display_list <- append(display_list, list(items))
    y <- y + nrow(items)
  }
  y <- y + 0.2
}
display <- bind_rows(display_list)
df_plot  <- display %>% dplyr::filter(!is_header)
headers  <- display %>% dplyr::filter(is_header) %>% dplyr::arrange(y)

# Colors
# 2) Paleta sincronizada con los nombres EXACTOS
section_cols <- c(
  "Totales"                           = "#A17C6C80",
  "Por causas (aj. por sexo y edad)"  = "#906F0080",
  "Edad (aj. por sexo)"               = "#5FA9B380",
  "Sexo (aj. por grupo de edad)"      = "#D2A899"
)

# Calculate limits, start from 0.5
min_lo <- min(df_plot$lo, na.rm = TRUE)
max_hi <- max(df_plot$hi, na.rm = TRUE)
xmin <- 0.5  # Force start at 0.5
xmax <- exp(log(max_hi) + 0.0 * (log(max_hi) - log(min_lo)))

# Rects
rects <- headers %>%
  transmute(
    section,
    ymin = y + 0.4,
    ymax = lead(y) - 0.4
  )
rects$ymax[is.na(rects$ymax)] <- max(display$y, na.rm = TRUE) + 0.6
rects <- rects %>%
  mutate(xmin = xmin * 1.005,
         xmax = xmax / 1.005)

# Labels position (moved left)
x_left_lab <- xmin * 0.5  # Moved further left

cap_h <- 0.25

p2 <- ggplot() +
  geom_rect(data = rects,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            inherit.aes = FALSE, fill = "#f7f7f7", color = NA, alpha = 0.6) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "#444444", linewidth = 1) +

  geom_segment(data = df_plot,
               aes(x = lo, xend = hi, y = y, yend = y, color = section),
               linewidth = 0.9) +
  geom_segment(data = df_plot,
               aes(x = lo, xend = lo, y = y - cap_h, yend = y + cap_h, color = section),
               linewidth = 0.9) +
  geom_segment(data = df_plot,
               aes(x = hi, xend = hi, y = y - cap_h, yend = y + cap_h, color = section),
               linewidth = 0.9) +

  geom_point(data = df_plot,
             aes(x = rme, y = y, fill = section, color = section),
             shape = 22, size=3,stroke = 0.6) +

  # Label over point with RME to 1 decimal, larger size
  geom_text(data = df_plot,
            aes(x = rme, y = y + 0.3, label = sprintf("%.1f", rme)),
            size = 5, color = "black", hjust = 0.5, vjust = 0) +

  geom_text(data = df_plot,
            aes(x = x_left_lab-.25, y = y, label = paste0("  ", label)),
            hjust = 0, vjust = 0.5, size = 3.6, color = "#222222", inherit.aes = FALSE) +

  geom_text(data = left_join(headers, rects[, c("section", "ymax")], by = "section") %>% mutate(y = ymax + .4),
            aes(x = x_left_lab-.25, y = y, label = section),
            hjust = 0, vjust = 0.5, fontface = "bold", size = 4.0, inherit.aes = FALSE) +

  scale_color_manual(values = section_cols, name = "Sección") +
  scale_fill_manual(values = section_cols, guide = "none") +
  #scale_size(range = c(2.6, 6.2), guide = "none") +

  scale_x_continuous(
    trans = scales::log_trans(base = exp(1)),
    breaks = c(0.5, 1, 2, 4, 8, 16),
    labels = function(x) format(x, big.mark = ".", decimal.mark = ",", trim = TRUE),
    expand = ggplot2::expansion(mult = c(0.3, 0.0))  # No expansion
  ) +
  coord_cartesian(xlim = c(0.5, xmax)) +  # Force xlim from 0.5
  labs(
    title = NULL,
    subtitle = NULL,
    x = "RME (IC 95%), escala logarítmica",
    y = NULL
  ) +
  theme_minimal(base_size = 15) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    plot.title.position = "plot",
    legend.position = "none"
  )

print(p2)

Warning in scale_x_continuous(trans = scales::log_trans(base = exp(1)), : log-2.718282 transformation introduced infinite values. log-2.718282 transformation introduced infinite values.

Code
if(grepl("ubuntu",pak::system_r_platform())){
  ggsave(paste0("forestplot_RME_log_corr.png"), dpi = 600)
} else{
  ggsave(paste0("forestplot_RME_log_corr.png"), dpi = 600)
}

Saving 7 x 5 in image

Warning in scale_x_continuous(trans = scales::log_trans(base = exp(1)), : log-2.718282 transformation introduced infinite values. log-2.718282 transformation introduced infinite values.

Session info

Code
#|echo: true
#|error: true
#|message: true
#|paged.print: true
message(paste0("R library: ", Sys.getenv("R_LIBS_USER")))

R library: /home/andres/Escritorio/SER2025/renv/library/linux-ubuntu-noble/R-4.4/x86_64-pc-linux-gnu

Code
message(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))

Date: 2025-10-14 00:08:33.054993

Code
message(paste0("Editor context: ", path))

Error: objeto ‘path’ no encontrado

Code
cat("quarto version: "); quarto::quarto_version()
quarto version: 
[1] '1.7.32'
Code
sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)|> 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Package' = 2, 'Version'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('R packages')),
      options=list(
initComplete = htmlwidgets::JS(
        "function(settings, json) {",
        "$(this.api().tables().body()).css({
            'font-family': 'Helvetica Neue',
            'font-size': '70%', 
            'code-inline-font-size': '15%', 
            'white-space': 'nowrap',
            'line-height': '0.75em',
            'min-height': '0.5em'
            });",
        "}")))
Code
#|echo: true
#|error: true
#|message: true
#|paged.print: true
#|class-output: center-table

reticulate::py_list_packages()|> 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Package' = 2, 'Version'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('Python packages')),
      options=list(
initComplete = htmlwidgets::JS(
        "function(settings, json) {",
        "$(this.api().tables().body()).css({
            'font-family': 'Helvetica Neue',
            'font-size': '70%', 
            'code-inline-font-size': '15%', 
            'white-space': 'nowrap',
            'line-height': '0.75em',
            'min-height': '0.5em'
            });",
        "}"))) 

Save

Code
if(grepl("ubuntu",pak::system_r_platform())){
  if (rstudioapi::isAvailable()) {
    # Code that needs RStudio
    wdpath<- dirname(dirname(rstudioapi::documentPath()))
  } else {
    # Code for non-RStudio environments (e.g., command line, server)
    # This part should use a portable method like knitr::current_input() or getwd()
    wdpath <- dirname(resolve_doc_dir())
  }
}
save.image(paste0(wdpath,"/congresosp2.RData"))
cat(paste0("Espacio de trabajo guardado en ",wdpath,"/congresosp2.RData\n"))

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
cat("Copy renv lock into cons folder\n")

if (Sys.getenv("RSTUDIO_SESSION_TYPE") == "server" || file.exists("/.dockerenv")) {
  message("Running on RStudio Server or inside Docker. Folder copy skipped.")
} else {
  source_folder <- 
  destination_folder <- paste0(wdpath,"cons/renv")
  # Copy the folder recursively
    file.copy(paste0(wdpath,"renv.lock"), paste0(wdpath,"cons/renv.lock"), overwrite = TRUE)
  message("Renv lock copy performed.")
}

Renv lock copy performed.

Code
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
time_after_dedup2<-Sys.time()

paste0("Time in markdown: ");time_after_dedup2-time_before_dedup2
Espacio de trabajo guardado en /home/andres/Escritorio/SER2025/congresosp2.RData
Copy renv lock into cons folder
[1] "Time in markdown: "
Time difference of 16.21523 days